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I. INTRODUCTION 

Much interest has recently been directed towards under- 
standing many-body dynamics in open systems away from 
thermal equiUbrium. This subject is not new, as the analogies 
between threshold phenomena in dynamical systems, such as 
the laser, and the conventional phase transitions have been rec- 
ognized over 40 years ago. However, recent experiments with 
ultracold atoms in optical cavities offer intriguing possibili- 
ties to explore the physics of strongly interacting atom-photon 
systems far away from thermal equilibrium from a new van- 
tage point. Many fundamental concepts of condensed matter 
physics, ranging from classification of phase transitions to the 
universal behavior of correlation functions in the vicinity of 
quantum critical points in the presence of driving and dissipa- 
tion, need to be revisited in Ught of these developments. 

In this paper, we investigate non-equilibrium phase tran- 
sitions for driven atomic ensembles interacting with a cavity 
mode that is subject to dissipation, focusing specifically on 
the dynamical superradiance transitions and associated self- 
organization of the atoms observed in Refs. 1-3. Due to the 
interplay of external driving, Hamiltonian dynamics and dissi- 
pative processes, the observed Dicke superradiance transitions 
[1,2] exhibit several properties [4—7] which are not present 
in the closed Dicke model [8-11]. Recently, it was shown 
that other more interesting quantum many-body phases, such 
as quantum spin and charge glasses with long-range, ran- 
dom interactions [12-14] mediated by multiple photon modes, 
could potentially be simulated with many-body cavity QED. 
Our longer-term objective is to develop a complete theory of 
these Markovian quantum glasses, including the dissipative 
processes and glassy phenomenology specific to quantum op- 
tics [15]. The present paper establishes the foundation for this 
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(focussing on the simpler physics with a single cavity mode) 
en route making a few smaller discoveries that have not been 
expressed clearly in the literature to date. 

In what follows we first review the physics of the non- 
equilibrium Dicke transition in optical cavities with con- 
ventional techniques of quantum optics. Using Unearized 
Heisenberg-Langevin equations, we demonstrate that at low 
frequencies, close to the phase transition, the system evolves 
into a thermal state with a high effective temperature, pro- 
portional to the atom-photon interaction strength. The static 
and dynamic critical exponents at the high-temperature phase 
transitions are analogous to those in a conventional laser (or, 
more precisely, optical parametric oscillator) threshold. 

To treat the interplay between extemal driving, dissipa- 
tion and many-body interactions in a more general setting, 
we next develop a unified approach to describe phase tran- 
sitions in open quantum systems based on the Keldysh path 
integral [16-22]. This approach is used to analyze both 
the low-frequency and the high frequency dynamics of the 
driven Dicke model with cavity photon loss and atomic spon- 
taneous emission. We re-derive the thermal properties of low- 
temperature dynamics, while at higher frequencies, we show 
that the photon and atom distribution functions exhibit quan- 
tum and non-equilibrium regimes. We beUeve the Keldysh 
calculations are not more involved, and sometimes simpler, 
than those of the usual quantum optics frameworks [23, 24]. 
At the same time, they faciUtate an easy comparison to other 
phase transitions of condensed matter physics. One of the 
objectives of the present paper is to make the Keldysh ap- 
proach more accessible to the broader quantum optics com- 
munity. At the same time, we hope that the Keldysh perspec- 
tive will be helpful for condensed matter physicists to un- 
derstand driven dissipative atom-photon systems-especially 
in view of the quahtatively different energy scales and bath 
properties in quantum optics. 

We expect the Keldysh approach to be directly applicable to 
other dissipative models such as the recently discussed central 
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spin model [25] or fermionic lattice models [26-29]. 

The paper is organized as follows. In Sec. 11 we introduce 
the Dicke model and perform a brief analysis with linearized 
Heisenberg-Langevin equations pointing out that the relevant 
low-frequency correlations are thermal. 

In Sec. Ill, we map the operators of the master equation 
and the associated LiouviUians for dissipative processes to the 
field content of an equivalent real-time, dissipative Keldysh 
action S{4>,a*,a] with (/> the atomic fields and a*, a the pho- 
ton field variables. As a warm-up, we first treat the case of a 
single, decaying cavity mode. Using this simple example, we 
present in Subsections III A - III D the program that we later 
in the paper repeat for the photon sector of the Dicke model 
with cavity photon loss (Sec. IV), and finally also for the atom 
sector including a simple model for atom decay (Sec. VI). 



II. THERMAL NATURE OF OPEN DICKE TRANSITION 

The Dicke Hamiltonian [8, 9] describes two-level sys- 
tems or "qubits" represented by PauU matrix operators cr?, cr^ 
coupled to a quantized photon mode represented by bosonic 
creation and annihilation operators , a: 

H = ojoa^a+'^j]cr^, + ^j]cr^{aUa) . (1) 

!=1 1=1 

Here coq is the photon frequency, co^ the level- sphtting of the 
qubits, and g the qubit-photon coupling, assiraied to couple 
all qubits uniformly to the photon. Eq. (1) is invariant under 
an Ising-type Z2-transformation, a — > -a and erf -erf. In 
the thermodynamic limit — > oo, and for sufficiently strong 
qubit-photon coupling g, the ground state of Eq. (1) sponta- 
neously breaks this Ising symmetry and exhibits a phase tran- 
sition to a "superradiant" phase with a photon condensate (a). 

In the context of ultracold dilute gases in optical cavities, 
Dimer et al. [4] proposed to implement the qubits using two 
hyperfine states of the atoms and showed that, close to the 
transition, the relevant Hilbert space can be exactly mapped 
to Eq. (1). Inspired by the work of Dimer et al. [4], the 
qubit states of the Dicke model were reaUzed using two collec- 
tive motional degrees of the Bose-Einstein condensate of the 
atoms in the cavity [2, 30], see [31] for a review. In that case 

becomes a collective recoil frequency and the two Dicke 
states are components of a dynamically forming charge den- 
sity wave. This open realization of the Dicke model in optical 
cavities with pumped atoms is different from the closed sys- 
tem Dicke model Eq. (1) for three principal reasons, which we 
now explain one-by-one: 

1. Coherent drive: The mapping to the Dicke model is per- 
formed in a frame that rotates with the (optical) pump fre- 
quency of the driving laser. In this rotating frame, the ex- 
plicit time-dependence of the original Jaynes-Cummings type 
Hamiltonian [54] is "gauged" away. For the effective parame- 
ters of the Dicke model Eq. (1), this leads to a shift of the cav- 
ity and atomic frequency scales: in particular ojq = iOc- cOp 
with cl>c the bare cavity frequency and aip the pump frequency 
(cf. Appendix A). This "hidden" time-dependence of the ef- 



fective parameters of the Dicke model turns the electromag- 
netic vacuum surrounding the cavity into a Markovian non- 
equilibrium bath. 

2. Cavity photon loss: In addition to the coherent dynam- 
ics generated by the Hamiltonian Eq. (1), there is a dissipative 
contribution consisting of cavity loss and spontaneous emis- 
sion for the atoms. The complete dynamics of the system den- 
sity matrix is modeled by a Markovian master equation, 

d,p = -i[H,p] + 2 '<a{2LapLl - {L%,p}) , (2) 

a 

where the curly brackets {, } denote the anti-commutator and 
La is a set of Lindblad or quantum jump operators, where a la- 
bels different degrees of freedom as well as different positions. 
The second piece of (2) describes the dissipative dynamics in 
terms of a Liouville operator X in Lindblad form, acting on 
the density matrix p. Here we have two contributions: cavity 
photon loss and atomic spontaneous emission. The former is 
modeled by the Liouvillian: 

-Ccav P = K{2ap^ - {a^a,p}) , (3) 

where /f is an effective decay rate (inverse Ufetime) of a cavity 
photon of the order a few MHz [2] . 

3. Atomic spontaneous emission: Modelling the dissipative 
dynamics of the atoms depends on the specific implementa- 
tion of the driven Dicke model usually involving local pro- 
cesses of each two-level atom separately. In Sect. VI, we 
account for spontaneous emission in an approximate way by 
resorting to a simphfied effective low frequency model. 



A. Heisenberg-Langevin analysis 

In this section, we study the effective dynamics of the Dicke 
model near the phase transition using conventional quantum 
optical techniques, namely, the Heisenberg equations of mo- 
tion for the atomic and photon operators. We begin by ne- 
glecting spontaneous emission of the atoms, but include cav- 
ity decay into a Markovian reservoir. In this case the dynamics 
is described by a set of operator equations: 

a = -icona-Ko, i=crf + T, (4) 

Vat 

a-^ = ico,crt - -^a^ia + a^), (5) 
cr} = -2^(crl-crT)(a + a'^'). (6) 

Va^ 

Here the force ;F - Tit) is a stochastic operator satisfying 
{TitWHt')) = 2K6{t - f) and {r{t)'^r{t)) = 0. This term 
is needed in order to preserve the commutation relation 
[a{t),aHt)'\ - 1, which would otherwise exponentially 
decrease. To analyze the dynamics below superradiance 
threshold, we assume that the atoms are fully polarized 
- |2,cr^ ~ -Nil, and neglect non-linear terms in 
the equations of motion. The resulting operator equa- 
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tions can be solved exactly in the Fourier domain as has 
been done in detail in the work of Dimer et al. [4] and 
we will not repeat their calculations here. We just add 
one simple point to their comprehensive analysis, namely 
that the relevant low-frequency dynamics of the photons 
occurs in the presence of a finite effective temperature. 
This is most easily seen from the equation of motion for 
the real-valued phase-space coordinate x{o}) = ia{oj) + 
aH-oj))/ V2wo- Defining stochastic force operators f(oj) = 

[T(cli) (k - ;(wo + cj)) + TH-io) (k + i^coo - oj))], 
can write the equation of motion as 



we 



i-K + ioif + a>Q ■ 



x{co) = f{io) . 



(7) 



The force operator satisfies 



Uf{oj)f{oj') + fW)f{o})) = K ^ 6{u + cj') . (8) 



(JJQ 



At low frequencies,we can neglect high-order terms in w. 
Eq. (7) becomes identical to the Langevin equation of a classi- 
cal particle in a harmonic potential with oscillation frequency 
a, defined by 



a - K + iOr,- 



OJ, 



(9) 



and friction constant 2k. This oscillator becomes dynamically 
unstable at the Dicke transition, where the frequency a van- 
ishes, given by the critical coupling 



gc 



4wo 



(10) 



The correlation function of the stochastic force operators on 
the right-hand- side of Eq. (8) becomes identical to the "noise 
correlations" provided by an equilibrium bath at a finite tem- 
perature 



4wo 



CJ, 



(11) 



However, in contrast to the temperature in an equilibrium 
problem, the low-frequency efl^ective temperature here is not 
a global property of the system, but is in general observable 
dependent. In Sec. IV D, we present a systematic, gener- 
aUzable way to extract low-frequency effective temperatures 
based on observable-dependent fluctuation-dissipation rela- 
tions. We already here quote the effective temperature for the 
atoms (see Sec. VI for the details of the computation): 



nefF 



■0)t 



4o), 



(12) 



where oj^ is the recoil energy and y an effective single atom de- 
cay rate. Because T^^ + Tf , the different parts of the driven 
system thus do not equilibrate to each other and, although the 



dominant low-frequency correlations are thermal, the system 
is not in a global thermal state. We remark that our definition 
of effective temperature does not coincide with those com- 
monly used in laser theory, as discussed in Sec. IV D: the 
former relates the fluctuations of the field to its response in 
the rotating frame, while the latter compares the fluctuations 
of the field to an equilibrium situation in the lab frame. 

We now demonstrate that the key critical exponent charac- 
terizing the non-equilibrium Dicke transition is dominated by 
the thermal modes in the spectrum. To that end, we compute 
the "photon flux exponent" [5] for photon occupation in the 
vicinity of the transition. The photon occupation is 



2(n) -I- 1 = {a^a + aa^) = 2ojo{{x^) + (/)) 



2ojf] 



1 + 



(13) 



where the expression for p = i{a - a^)/ ^J2ojq can be obtained 
in complete analogy to that of x. We can compute (x^) using 
the equilibrium partition function [32, 33] Z = exp(-F/T^^^), 
with F = ^a^x^. Performing the Gaussian integral yields 



2(n)-i-l =21 + 



,.2 \ 



rtf[ 



1 



ojo- 



(14) 



with the static photon flux exponent Vx = 1, as already noted 
in Ref. 5, 6. There the authors also noted the discrepancy 
with the zero temperature case [11], where (n) ~ I /a ~ 
l/\g^ - gl\^^^ and = 1/2. Here we offer an explanation 
of this discrepancy: the dominant photon modes at the Dicke 
transition have a non-zero effective temperature. 

The above results indicate that, from the point of view of 
phase transitions, it is incorrect to call the driven Dicke tran- 
sition a quantum phase transition even though it is "made of 
quantum ingredients" (two collective motional states of a BEC 
[2]), and the photon spectrum contains a squeezed quantum 
mode. Instead, it should be regarded as a classical phase tran- 
sition governed by the properties of the dynamical universal- 
ity class of the classical Ising model, the so-called "Model 
A" of Hohenberg and Halperin [34]. The effect that dissi- 
pation induces a finite effective temperature is not new. In 
several other condensed matter systems [35-37], the coupling 
to a non-equilibrium bath typically admixes the pure many- 
body states of the closed system, transforming pure quan- 
tum phase transitions into thermal phase transitions (see also 
[17, 18, 25, 38, 39]). What is perhaps more surprising is that 
neither the low-frequency effective temperature for the pho- 
tons nor for the atoms is set by the cavity loss rate k, but in- 
stead set by the atom-photon interaction g and the effective 
single atom parameters, respectively. 

Being complex-valued objects, the photons have two nor- 
mal modes. As we saw, one quadrature is thermally amplified 
and diverges at the Dicke transition. Its orthogonal quadrature 
is quantum squeezed [4], remains gapped at the transition and 
therefore does not influence the thermal nature of the phase 
transition. These attenuated and amplified quadratures arise 
naturally as the eigenmodes of the photon correlation func- 
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tion (see Sec. IV B). 

The remainder of the paper is dedicated to the development 
of a unified, and generalizable, Keldysh approach. The Dicke 
model will be used as a prototypical test object and we com- 
pare oiu' results to those of other approaches, where available. 



in. KELDYSH APPROACH FOR CAVITY VACUUM 

In this section, we introduce the real-time Keldysh formal- 
ism and fix our notation by considering the case of a single 
electromagnetic mode in an open cavity (without atoms). We 
also indicate how to relate the operators of the master equa- 
tion (2) to the field content and choice of time contour in a 
Keldysh action (see also Refs. [16, 19, 40]). 

The decay of a single boson (the cavity photon) into a con- 
tinuum of modes (the external vacuum) is described by the 
master equation 

d,p = -i[ojoa^d,p] + K(2apa^ - {a^d,p}). (15) 

This equation results from the (fully unitary) Heisenberg 
equation for the coupled system-bath setting, where the sys- 
tem is described by the degrees of freedom a, d^ and the bath 
by a continuum of harmonic modes. Eq. (15) is obtained by 
eliminating ("integrating-out") the bath in the Born-Markov 
and rotating wave approximation (cf. e.g. [23]). The integra- 
tion of the bath variables gives rise to an effective evolution 
including dissipative terms. Performing the same program in 
path integral formulation, we obtain the Markovian dissipa- 
tive action 

dt{a\{idt - o}(i)a+ - aliid, - a)o)a- (16) 

oo 

-iK[2a+a*_ - {a\a+ + a* fl_)|j. 

In the path integral formalism, the quantum mechanical opera- 
tors are replaced by fluctuating, time-dependent and complex- 
valued fields (we omit the time argument for notational sim- 
plicity). The fact that the density matrix can be acted on from 
both sides, as reflected in the Heisenberg commutator struc- 
ture of the original evolution equation, finds its counterpart in 
the presence of a forward (-I-) and backward (-) component of 
the fields. The former is associated to an action on the den- 
sity matrix from the left, and the latter to the right. Indeed, in 
the first of line of Eq. (16) there is a relative minus sign be- 
tween the terms involving the two components, reflecting the 
Heisenberg commutator structure of Eq. (2). The terms in the 
second line instead display the characteristic Lindblad form; 
the "jump" or "recycling" term is represented by an explicit 
coupling of the two contours. 

It is convenient [20-22] to introduce "center-of-mass" and 
"relative" field coordinates, - (a+ + «_)/ V2, - («+ - 
a- )/ V2. These new coordinates are often referred to as "clas- 
sical" and "quantum" fields, because the first can acquire an 
expectation value while the second one cannot in the absence 
of sources. In this basis, and going to frequency space, we 



write 

where we used the notation = ^, and adqit) = 
e^'"'aciq{(jS). This classical-quantum basis is often referred 
to as RAK basis: the entries are the inverse Retarded (lower 
left) and Advanced (upper right) Green's functions, and the 
inverse Keldysh component. The RAK action Eq. (17) can be 
easily inverted to deliver the photonic Green's functions 

(18) 

where the Keldysh Green's functions is a matrix product 

G'^ioj) = -G\aj)D'^{oj)G\aj) . (19) 

For the open cavity of Eq. (16) the RAK inverse Green's func- 
tions are 

[G-'f'^=oj-oj,+I.ll\ i)^=S^ (20) 

with the "self-energies" 

S^ = -j*:, 'Ll = +iK, l^ = 2iK. (21) 

It is a key property of a Markovian system that the Keldysh 
component in Eq. (21) is frequency independent. As can 
be seen from App. A, this is due to a separation of scales be- 
tween (i) the large pump (ajp) and cavity frequency (w^), both 
optical frequencies in the Tera Hertz range (~ lO'"* Hz cor- 
responding to temperatures T ~ 10'^ K), and (ii) the charac- 
teristic frequencies of the electromagnetic vacuum outside the 
cavity (< 10'^ Hz corresponding to temperatures T < 300K). 

In the literatiu'e, it is often argued that a frequency inde- 
pendent, nonzero inverse Keldysh component indicates an ef- 
fective finite temperature state. The Markovian lossy cavity 
is a simple counterexample: Even though the inverse Keldysh 
component is constant ~ 2k, the state is pure and the effective 
temperature zero, as we will argue below. 

We next introduce the key propagators that encode the sys- 
tems' response and correlation functions. In equilibrium, the 
two are rigidly related by the Bose (or Fermi) distribution 
function; out-of-equilibrium, no such a priori knowledge is 
available, and it is important to distinguish them. 



A. Cavity spectral response function 

The spectral response function encodes the system's 
response to active, external perturbations such as time- 
modulated extemal fields coupling to spin or charge operators, 
for example. The spectral response function is the difi'erence 
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between the retarded and advanced Green's function: 

J?l(w) = i(G*(w) - G\a>)). (22) 
In the scalar case considered here, we have 

Ji,A<^) = -llmG^iu), (23) 

The frequency-integrated spectral response function is nor- 
malized to unity, because of the exact commutator relation of 
the bosonic degrees of freedom 



do) + 



(24) 



This "sum rule" is an exact property of the theory valid in- 
and out-of- equilibrium. In our example of one cavity mode. 



2k 



(25) 



respect to the frequency ^^^^(-a)) = -^^eq^^)' while in the 
Markovian system (21) it is even. 
Using Eq. (28), we obtain 



J?l^'^(w) = hm 



2e 



aa* ^ ^ 



2e 



hm 



x9 T coth ^ 

(w - Wo)2 + e2 2T 



= 2nd{cj - wo), (29) 



In this noninteracting case, the spectral response function is 
fully centered at the isolated mode with frequency ojq. We 
observe that formally, the thermodynamic equilibrium limit 
can be seen as a situation with an infinitesimal loss, replacing 
K ^ e, and the replacement in the inverse Keldysh component 
2i/c^2iecoth^. 



D. Cavity distribution function and low-frequency effective 
temperature 



B. Cavity correlation function 

The correlation function encodes the system's intemal cor- 
relations, for example the frequency-resolved photon spec- 
trum of the intracavity photon fields. In the steady state, the 
photon correlation function is related to the Keldysh Green's 
function by 

Caa<t) - mt\a\0)}) = mra\0) + a\0)a{t)) - iG^(f). 

(26) 

Here the last identity is valid only in the specific case of a 
scalar Keldysh Green's function. At equal times this relation 
results in 

/J 
-^G\a}). (27) 

For the single decaying cavity mode, characterized by 
Eqs. (17-21), it is easy to show that Ca„;(w) - y\„„-<{(jS) and 
the frequency integral over the Keldysh Green's function is 
unity yielding {a^a) = 0. As expected, the steady state corre- 
sponds to the cavity vacuum. 



C. Comparison vnth a closed system at equilibrium 

In the absence of dissipation, the cavity becomes an iso- 
lated harmonic oscillator. Its inverse Green's functions are 
still given by (17) with the self-energies serving only as regu- 
larization parameters. At equilibrium. 



= —le. 



S^eqC'^) - 2iecoth 



(28) 



Here e — > at the end of the calculation and T is the actual 
temperature. Note that the Keldysh component is odd with 



The response and correlations allow us to define a 
fluctuation-dissipation relation, by introducing the distribu- 
tion function F{oj): 



G^icS) = G\o))F{io) - F{io)G\o}) 



(30) 



{G\cS)\'f{cS) - F{cS){G\oj)\\ 



where the equivalence holds due to Eq. (19). At thermal equi- 
librium the distribution F is universal and equals to the unit 
matrix times 



Fj%{oj) 



!7EQ 



coth^ =2n2,(fJ) + l, 

sign(w), 
2T 

0} 



+ ... 



(31) 



with the Bose distribution riBix) - (exp x - 1) ^ The unit ma- 
trix in field space signals detailed balance between all subparts 
of the system. 

In the present case, the system is out-of-equilibrium due to 
its driven and dissipative nature. A notion of a temperature 
is not a priori meaningful: Neither must the driven system 
equilibrate to an external heat bath with temperature T (in our 
case, due to the separation of scales underlying the Markov 
approximation, this temperature would be effectively 7 = 
compared to the system scales), nor do the different subparts 
of the system have to equilibrate with respect to each other. In 
this work, we argue that a notion of a temperature nevertheless 
emerges as a universal feature of the low frequency domain 
of Markovian systems. It is introduced by computing the F 
matrix through Eq. (30) and comparing the low-frequency be- 
havior of its eigenvalues with the equihbrium result of Eq. 
(31). In particular, if F has a thermal infrared enhancement 
~ l/(jj for small frequencies, its dimensionful coefficient is 
identified as an effective temperature. This notion of a "low 
frequency effective temperature" (LET) becomes particularly 
relevant in the vicinity of a phase transition, where the spec- 
tral weight encoded in G*, G^ is concentrated near zero fre- 



6 



quency. Below, we will use this concept to establish a con- 
nection of Markovian quantum systems to the classical the- 
ory of dynamical universaUty classes according to Hohenberg 
and Halperin [34]. Moreover, we find that, while all governed 
by the 1 /a» divergence in the distribution function, diff'erent 
subparts of the system exhibit different LETs. In contrast to 
the global temperature present in thermodynamic equilibrium, 
the LET is not an external parameter but rather a system im- 
manent quantity, determined by the interplay of unitary and 
dissipative dynamics. 

In case of a decaying cavity, the Green's functions are 
scalars and we can easily invert (30) to obtain 



1. 



(32) 



confirming that the cavity vacuum has a zero effective tem- 
perature. Moreover, as for a pure quantum state in the equilib- 
rium case at r = 0, the distribution function here also squares 
to a unit matrix, F^(ai) - 1. 

An important difference between the zero temperature equi- 
librium case and Markovian case appears in the sign of the 
distribution function. In the former case (as for any equilib- 
rium distribution), F{a)) is anti-symmetric with respect to the 
frequency at. On the contrary, for a Markovian bath the distri- 
bution function is symmetric with respect to (o; one signature 
of a strongly-out-of-equihbrium system. 



IV. KELDYSH APPROACH FOR PHOTON OBSERVABLES 

We now analyze the Dicke model Eq. (1) with the path 
integral approach explained in the previous section. We in- 
clude cavity photon loss but defer the inclusion of sponta- 
neous emission to Sec. VI. Assuming homogeneous qubit- 
cavity coupling, one can use the large-A^ strategy of Emary 
and Brandes [11]. We introduce collective large-A^ spin op- 
erators = i li^li cr^ and 5^ = i + o",") ' to write 
the Dicke model (1) in terms of one large spin coupled to the 
cavity photon mode. 



H 



' (a"'" + a) 



(33) 



We then express the spin in terms of a Holstein-Primakoff 
[11, 41] boson operator b, defined by 5j. = -N/l + b^b, 
5+= yfW^P « V^(l-n/(2Af))^"',and5* = (5++5-)/2. 
Neglecting unimportant constants we obtain the normal or- 
dered Hamiltonian: 



H = a>Da^a + cj,Pb 

+g{a + a^){b + P - ^b\P 



(34) 



At A' — > oo, the last, non-quadratic term vanishes and the 
problem reduces to a hnear system of two coupled bosonic de- 
grees of freedom one of which (the cavity mode a) decays into 
a Markovian bath. As outlined in section III, we can trans- 
form the Liouvillian Eq. (3), with Hamiltonian Eq. (34), into 



an equivalent Keldysh action with 

S = Sa+Sb+Sab 



(35) 





CJ- CJ, 







bd 



where 5 a is given by Eq. (17) and the interaction in terms of 
the "classical" and "quantum" fields reads 

Sab = -5 r(l - 2k)[(«e + + Ki) + («c/ + aa){bq + fop] 

■Jo) 

- ^{ [iaci + KiKb, + fop + (a, + ap(fo,, + foj,)] [fo>,, + fo;fo,] 

+ [ia„ + al,)ib„ + bl,) + (a, + apib, + fopj [fo^^fo, + fo;forfl | . 

We now demonstrate that the static saddle point solutions of 
this action reproduce the results of other approaches [4]. Vary- 
ing 5 with respect to the quantum components of the fields 
and substituting adit) - ao, bdit) = foo, a« = 0, bq - 0, we 
obtain the coupled equations 



dS 



1 



(-(Oo + iK)ao -g^l--b'o\2bo = 0, (36) 



dS_ 
5fo* 



-Wjfoo - ^ ( 1 - -foo I (ao + ap = , 



(37) 



where we chose foo = b*y As expected, these saddle-point 
equations admit solutions with non-zero "ferromagnetic" mo- 
ment foo, and super-radiant photon condensate ao, 



fo, 



'0 




ao 



g"- OJo - IK 

for atom-photon couplings larger than a critical value 



(38) 



8> 8c = 



4a>o 



(39) 



We are now going to integrate-out the atomic field fo and ob- 
tain an effective action describing the photons in the normal 
phase (g < gc), where ao = foo = 0. In App. B, we give the 
corresponding expressions in the superradiant phase. 

In the limit of N ^ 00, we can safely neglect the terms 
proportional to 1 /N in (35). We can re-write Eq. (35) as 8 X 8 
matrix multiplying the 8-component fields: 



Odioj) 
bd(co) 

a*q(-a)) 

bqiCO) 

b'qi-co) ) 



(40) 
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with the action 



(41) 



The subscript N stands for normal phase and the dagger t de- 
notes transposition and complex conjugation. The block en- 
tries are 4x4 Green's functions given by 



O) - ajd + iK -g -g 

-W - Wo - IK -g -g 

-g -g a)-cD, 

-g -g -w-w, 

Z)f^4 = 2i diag(/f,/f,0,0). (42) 

To obtain the photon-only action, we now integrate out 
the Holstein-Primakoff field b and get a Keldysh func- 
tional integral that goes only over the photon fields: - 
jD{a*, aje'^^p'™'"*'"^. with the photon-only action 



•5 photon [t^ ,0]— I ^/ 



(W) 



A4(lo) . 
(43) 



The photon four-vector collects the classical and quantum 
field components 



A4(w) 



flrf((jj) ^ 
fl*,(-w) 

aq((D) 



(44) 



and the block entries are 2 x 2 photon Green's functions which 
we now analyze one-by-one. 



A. Photon spectral response function 

The inverse retarded Green's function of the photons is 



2x2- 

OJ - OJQ + IK + 2*((x») 

[s*(-w)]* -co - Wo - iK + [l*(-w)]* 



where the interaction induced photon-self energy reads 



S«(w) = . 



(dp- — (dp 



(45) 



(46) 



The characteristic frequencies of the system are defined by 
the zeros of the determinant [G2X2] ' ('^)' corresponding to the 
poles of the response function G^^^ibS). Due to the symmetry 



(47) 




FIG. 1: (Color online) Schematic plot of the position of the poles of 
the retarded Green's function Eq. (48). (a) At zero coupling g = 0, 
two poles can be associated with the photonic mode ai = ±010 + Ik 
and two with the atomic mode uj = ±(1)-. (b) In the presence of 
a finite coupling g, the modes hybridize and the corresponding fre- 
quencies are shifted in opposite directions, (c) When approaching the 
transition, two solutions become purely imaginary and correspond to 
damped modes, (d) At the transition point g = gc one of the poles 
approaches zero, making the system dynamically unstable. 



the poles come in pairs, such that {A} = {-/?*), meaning that 
they either are pure imaginary or come in pairs with opposite 
real part. The explicit solution of 







1 



Wo -H 



2g^co, 



-(to + iK) 



2g^oj, 



det[Gf^2(«)] 

(48) 

yields four poles, schematically plotted in Fig. 1. Note that, in 
the vicinity of the phase transition, two poles become purely 
imaginary. This phenomenon seems to apply to generic dissi- 
pative transitions as we further describe in App. C. Indeed it 
has been previously observed for a dissipative critical central 
spin model [25]. Overdamping of collective modes, due to a 
similar mechanism, has also been found in dissipative multi- 
mode systems in symmetry broken phases [17, 18]. 

The imaginary part of the first diagonal element of G2^2(^) 
corresponds to the photon spectral response function ^aot, 
defined in (23), and is plotted in Fig. 2 for different values 
of the coupling g. In the absence of atom-photon coupling 
g = (dotted curve) there is a single resonance peak at fre- 
quency (JL> - Wo, broadened by the cavity decay rate k. A finite 
coupling g (dashed curve) "collectively Rabi-splits" the reso- 
nance [42] in two distinct peaks, corresponding to two distinct 
poles of the system. Upon approaching the Dicke transition 
(solid curve), the spectral weight is shifted towards the low- 
frequency pole; a precursor to the superradiant cavity mode. 
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FIG. 2: (Color online) Photon spectral response function A^^t (li') and 
correlation function C„„t(w) as a function of real frequencies cj. Nu- 
merical parameters: ojq = w, = 1, = 0.2 (leading to ~ 0.51), 
and g=0 (dotted), 0.25 (dashed), 0.5 (solid). 



laser drive process compensating for the loss. 

The matrix structure of the Keldysh Green's function 
G\^2^(jS) can be conveniently exploited to compute the quadra- 
ture fluctuations for a general phase angle Q 



where xg is defined by 



(50) 



(51) 



The corresponding equal time, frequency-integrated fluctua- 
tions {xg(t)xg(t)) - {xg(cij)xg(-tij)) are plotted in Fig. 3 and 
diverge at the Dicke transition for all angles 6 except for 0* 
defined by 0* = tt - tan"'((Wo/A-). The angle 6* can also be 
obtained as the phase angle for the non-diverging eigenmode 
of the zero-frequency limit of the Keldysh correlation func- 
tion, thereby naturally yielding the attenuated and amplified 
quadratures alluded to in Sec. 11. 

Here in the case of the driven Dicke model, the equal-time 
photon fluctuations of xg* inside the cavity are independent of 
the atom-photon coupling and in particular are not attenuated 
below the vacuum noise level (the g - Q limit, without atoms 
in the cavity). This is different from the case of the optical 
parametric oscillator [43] where, at threshold, the equal-time 
fluctuations of the non-diverging intra-cavity quadrature are 
reduced to 50% of the vacuum level. This difference can be 
traced back the different frequency dependencies of the effec- 
tive driving term in both situations. In the parametric oscil- 
lator, the driving of the cavity occurs via a classically-treated 
photon pump laser, and the driving amplitude is typically set 
to a constant coherent field amplitude. In the present case of 
the driven Dicke model, the effective driving of the cavity is 



B. Photon correlation function 



The inverse Keldysh component of the action Eq. (43) is 



D 



2x2 



2iK 
2iK 



(49) 



Green's function Gf 2('^) 



t(w) is a 2 X 2 matrix. 



and the Keldysh 

-^2x2('^)^2x2^2x2 

The first diagonal element of Gf^^, Ca^t = i(l 0)G^^2(1 0)^ 
coiTesponds to the photon coiTelation function defined in 
Eq. (26) and is plotted in Fig. 2. As noted above in Eq. (27), 
its frequency integral gives the steady-state photonic occupa- 
tion and it can be shown to diverge at the Dicke transition 
according to Eq. (13). This result will be explicitly derived in 
Sec. V A using a low-frequency effective description of G^. 
The photon number diverges in steady state despite the fact 
that the system undergoes photon loss, and no explicit photon 
pumping occurs within the model. The reason is that the cou- 
pling constant g is an effective parameter, which in any con- 
crete physical realization microscopically involves a coherent 



V 




1.0 0.5 



FIG. 3: Frequency integrated equal-time correlations {(x/)(t)) > as 
function of the angle d, for different values of the coupling strength 
approaching the transition at ~ 0.51: g=0 (dotted), 0.25 (dashed), 
0.5 (solid). The fluctuations in the quadrature 6' = ;r-tan~'(6Jo/*') ~ 
1.768 are independent of the coupling strength. Numerical parame- 
ters: ojq = ai, = 1.0, K = 0.2 
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mediated by the atoms via virtual absorption and emission of 
photons from the pump laser into the cavity. The correspond- 
ing driving term ~ g^a)^j{a)\ — w^) is maximal for frequencies 
0) of the order of the atomic detuning cj- and vanishes for large 
frequencies. 

Nevertheless, also in the driven Dicke model, the experi- 
mentally relevant homodyne spectrum of '^he cavity 
output field shows noise reduction below the vacuum level in 
the 0* quadrature [4]. Following the standard "input-output 
theory" [43, 44], it is possible to show that the homodyne 
spectrum can be linked to the Keldysh response function via 



2iK 
-2iK 



G2v2(<^) - l2x2 



12x2, 



(52) 



By applying the transformation (50) 



2j^2' it is then possible to compute the fluctuations of 



where ||M|r = M"'"M, 
to G™' 

the output quadrature {xoM,e{<ii)xont,e{-(^)) ■ This quantity has 
been studied in detail in Ref. [4]: at the Dicke transition, 
the zero-frequency component of 9 - 6* tends to the maxi- 
mally attenuated value of -1/(4(jJo). It should be noted, how- 
ever, that the equal-time, frequency integrated fluctuations 
(-^out.e* (f )'^out,e* (0) are zero and therefore not attenuated below 
vacuum level. 
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FIG. 4: (Color online) Positive eigenvalue of the distribution func- 
tion F{(jj) for the same parameters as in Fig. 2. At low frequencies 
the distribution diverges as IT^ff/o), where the finite effective tem- 
perature is proportional to the photon-atom interaction, Tejf ~ g~. 



Eq. (43) we find: 



(54) 



where cr, and cTx are Pauli matrices. The F matrix is hermitian 
and traceless, so its two eigenvalues are real and opposite: 



C. Comparison with a closed system at equilibrium 

In the case of a closed system at equilibrium, the retarded 
Green's function is the same as Eq. (45), up to the replace- 
ment A- — » e, and letting e — > 0. The corresponding spectral 
response function ^„gi{aj) is similar to the one shown for the 
Markovian case, with the narrow peaks substituted by delta- 
functions at the resonant frequencies. The Keldysh compo- 
nent of the inverse Green's function reads, at zero tempera- 
ture: 



D 



2x2, EQ 



2/esign(w) 

-2/esign(w) 



(53) 



Again note the different symmetry under frequency reflection 
with respect to (49). 



D. Photon distribution function and low-frequency effective 
temperature 

We now show that the above mentioned difference between 
the Markovian case and the equilibrium case leads to the gen- 
eration of a "low frequency effective temperature" (LET) for 
the former. For this purpose, we calculate the distribution ma- 
trix F, defined in Eq. (30). To this end, recall that at thermal 
equilibrium, F - coth(w/2r)l: it exponentially approaches 
unity at high frequencies (|w| » T), and diverges as 2T joj at 
low frequencies (|w| <s: T). For our Markovian problem, using 



1 + 



1 



1 - {(jj/to^y 



(55) 



Fig. 4 shows the behavior of the positive eigenvalue in two 
points of the phase diagram, nearby and far away from the 
transition. In both cases, at low frequencies the eigenvalue 
diverges as l/co. Exploiting the analogy to the equilibrium 
case, we obtain the LET 



! eff 



(56) 



We find that Teff is not proportional to the decay rate k. The 
temperature is rather proportional to the effective interaction 
between the spin and the photon, g is the scale that leads 
to a competition of unitary and dissipative dynamics, and the 
LET is a measure of this: For g = 0, the steady state of the 
dissipative part of the dynamics (empty cavity) is an eigenstate 
of the Hamiltonian, while this is no longer the case for any 
finite g. In the limiting case of ^ — » 0, the LET goes to zero. 

A closer inspection of (55) reveals an important difference 
between the Markovian bath and thermal equilibrium, related 
to the presence of a second energy scale, a>,. If a>~ » T^ff, this 
energy scale does not affect the crossover between the quan- 
tum and classical regimes, which then proceeds monotonously 
similar to an equilibrium problem. If, on the other hand, 
<s; Teff, the quantum-classical crossover of the Marko- 
vian bath occurs in an unusual way, highlighting the non- 
equilibrium nature of the problem. Starting from a divergence 
at zero frequency, the distribution function (55) decreases as 
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Tetf/w, in analogy to an equilibrium system at finite temper- 
ature. Then, instead of monotonously decreasing towards the 
quantum regime where f ^ 1 , it exhibits a second divergence 
at a> = Wj. Since the spectral weight vanishes sufficiently fast 
in this regime, the correlation functions still remain finite; the 
pole in F accounts for a different scaling of correlations and 
spectral properties in this regime (cf. Fig. 2). At higher fre- 
quencies, it finally tends to one, following the non-equilibrium 
curve / K Teff/{aj - w j. The approach to the quantum regime 
/ Si 1 is polynomial, unhke the exponential approach in the 
equilibrium case. In App. D, we show that the thermal l/co 
divergence, leading to a finite LET, is generic for Markovian 
systems. 

We note that our definition of an effective temperature is 
not the one commonly used in the context of laser theory [24]. 
In this context, the effective temperature is used only to de- 
scribe the fluctuations of the photonic field, as compared to the 
equilibrium fluctuations in the lab frame. As a consequence, 
the divergence of the photon number at the phase transition 
is always associated to a diverging effective temperature. In 
contrast, our low-frequency effective temperature (LET) de- 
scribes the ratio between the fluctuations and the response of 
the system in the rotating frame. It is finite at the transition 
and, as we will see, allows us to map the Dicke transition to an 
existing dynamical universality class of equilibrium systems. 

V. PHOTON LANGEVIN EQUATION FROM KELDYSH 
ACTION AND FINITE-W CORRECTIONS 



and the determinant of Gf^a^'^-' given by Eq. (48). 

As well known, any quadratic Keldysh action is equiva- 
lent to a linear Langevin equation. Starting from a generic 
quadratic action (57) one introduces a Hubbard-Stratonovich 
"noise" field f(aj) to obtain: 

(59) 

(This action is equivalent to Eq. (57), as can be expUcitly 
shown by performing the Gaussian integral over f{oS), and us- 
ing G^icS) - G^{—oS) - see Ref. 21, 22 for more details.) In- 
side the Keldysh partition function = ^ D{xq; Xd, f] e'^'-f, 
the integration over Xq then takes the form of a delta function 
with argument 

[Gi,(w)]-ixrf(w) = f(uj). (60) 

The remaining last part of the action Eq. (59) involves only 
f(aj) and can be thought of as the statistical weight of a Gaus- 
sian random variable with correlations 

</(«)/(w')) = -^Di(<^)5(aj + w') . (61) 

Eq. (60) then becomes a stochastic equation of motion for the 
bosonic field x(w) = Xdica), identical to the Langevin equation 
(7) obtained in Sec. II A. 



We now rederive the simple scalar Langevin equation (7) as 
the stochastic equation of motion of the Keldysh action. We 
also determine how the photon correlations scale with system 
size when approaching the Dicke transition. 



A. Low-frequency effective Langevin equation for Af ^ oo 

In this subsection we derive a simple description of the 
photon-only action (43), focusing on the x - {a + a*)/ y/laiQ 
quadrature, by mapping its Keldysh action to a linear stochas- 
tic equation. The resulting Langevin equation coincides with 
the equation of motion of a classical particle in an harmonic 
confinement, coupled to an equilibrium bath at finite tem- 
perature consistent with the Heisenberg-Langevin analysis of 
Sec. II A. 

To obtain the Keldysh action of x we replace a = 
Vwo/2(x + ip) and a* = -\/wo/2(x - ip) inside Eq. (43) and 
integrate out the p-component. This procedure gives: 



Sxx= I [xcii-oj) Xgi-Oj)) 

Jul 



-1 \ 



where 



+ 1^ + 0? 



[GiSoj)] 



[GUco)]-' = 



Xcl(co) 
Xq(a>) 

(57) 



det[GR^2(w)]' 



B. Dynamic critical exponent 

In addition to the photon flux exponent, we identify a sec- 
ond indicator of criticality, the dynamical exponent. This ex- 
ponent governs the decay of the two-time correlations close 
to criticahty. Going back to the Langevin equation (60) and 
keeping the co^ terms we obtain: 

[Vo/ + liKco + a^) x(aj) = /(w), (62) 

0)1+ + o/ 

{fico)fico')) = 6{oj + CO') , (63) 

where we defined a dimensionless parameter y = 1 + 

4(jj()g- / oj^- . For simplicity we further approximate V « 1 and 
obtain the correlation function [55] 

mt),xm)=iGUt) (64) 



In 2wo[(2/fw)2 + (Q'2-w2)2] 



(58) 



= \m^ {a)l + + a^] cos ( Vm^A 

+ Vm^^Wo + /c:^ - or^)sin( Vm^?)] , 

where we defined = a^-K^. In the vicinity of the transition, 

for < a < K, the frequency m becomes purely imaginary and 
the oscillatory behavior in the above expression disappears. 
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As already mentioned, this is a generic feature of dissipative 
phase transitions (see App. C). For sufficiently large times and 
approaching the transition a — > (where only the closest pole 
to zero contributes), we then obtain 



The correlation time: 



2k 



1 



(65) 



(66) 



\gc - gV' 

is governed by a dynamical exponent v, = 1 . 

C. Finite size scaling 

We now include non-linearities for the photon dynamics 
arising from a finite number of atoms A^. The most relevant 
term in the Langevin equation is a frequency independent cu- 
bic term: 



fit). 



(67) 



This equation defines the dynamical critical theory of an Ising 
transition with no conserved quantities, the so-called "Model 
A", for n - \ degrees of freedom and d - Q dimensions 
[34]. The frequencies a and fi need to be determined from 
the microscopic theory. The former has already been derived 
in Eq. (9). To obtain on the other hand, we demand that the 
effective Langevin equation reproduces the same saddle-point 
as the original action Eq. (35). Using x - {a + a*)/ yJlcDQ and 
Eq. (38) we then obtain: 




V2wo(g^ 



2) 



^^ = 2 5-. (68) 



The parameter /? describes deviations from the mean-field so- 
lution, due to the finiteness of the system and, indeed, it van- 
ishes in the thermodynamic limit — > oo. 



With Eq. (13) and {x^) = 

1^2 ^2 



and the effective 



free energy F{x) - ^a^x^ + jff'x'^, we can write the photon 
occupation at the critical point as: 



2<«>, + 1 



2w(: 



1 + 




r(i/4) 



Va^ (69) 



Here r(3/4)/r(l/4) ^ 0.338 is the ratio of two Gamma func- 
tions. We find that the photon occupation at the critical point 
scales like A^'^^. This result indicates a further distinction be- 
tween the present model and the zero-temperature quantum 
phase transition [45, 46], for which (n) ~ A^'^^. 

To substantiate our findings we compare them with the nu- 
merical simulation of the Dicke model (1) with cavity loss (3) 
using the Monte Carlo Wave-Function (MCWF) method [47], 
as implemented in the open-source C-I-+QED library [48]. 
The numerical results in Fig. 5 compare well to the effec- 
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FIG. 5: (Color online) Photon occupation in the vicinity of the 
Dicke transition for A'=5 (red), 10 (green), 20 (blue), 40 (magenta). 
Crosses: result of the numerical simulation of the full Master equa- 
tion. Dashed line: mean field result (14), valid in the thermodynamic 
limit N CO. Solid lines: elfective Langevin equation, including 
leading 1/A' corrections. The effective model is valid only around 
the critical point and its validity region narrows with increasing A'. 
INSET: log-log plot of the photon occupation at the critical point. 
Crosses: numerical results. Solid line: effective model (69), predict- 
ing a A'''- dependence. Model parameters: lo- = 2.0, oiq = 1.0, 
K = 1.0, giving = 1.0. Other parameters used: number of trajecto- 
ries Miaj = 10, time step dt = 1, total time T = 400. 



live Langevin approach (69) and approach the predictions of 
the mean field approach (14) for A^ — > oo. We emphasize that 
no fitting parameters were used when comparing the numerics 
with the effective model. Independent calculations [49] indi- 
cate that the critical exponent of the photon occupation may 
differ from the predicted rational scaling dimension. A down- 
shift from our analytical result 0.5 to 0.41 + 0.02 is found. 
This value is compatible with our present numerical results as 
well. This discrepancy is probably due to fact that we extract 
the critical exponent at leading order in the 1/A^ expansion, 
which may be a hint at anomalous critical behavior. Such an 
anomalous dimension can occur at equilibrium as well (both 
at zero and finite temperature) and does not probe a priori 
the non-equilibrium nature of the transition. However, as the 
limit A^ — > oo leads to an approach of a Gaussian fixed point, 
while any finite value of A^ corresponds to a non-critical sys- 
tem, such a result is surprising. A more detailed investigation 
is left for the future. 

We now consider the dynamical properties (i.e. non-equal- 
time coiTelations) of the open Dicke model at the critical point 
and show that they differ from those of an equilibrium sys- 
tem at finite temperature [8-10]. To that end, we study the 
Langevin equation (67) at the critical point a = under the 
scaling transformation /3 — > (cf. Eq. (68)). How should 

we rescale f, x, and / such that the Langevin equation (67) is 
still satisfied? First we note that, due to the presence of a finite 
(effective) temperature (/(?)/(?')) - T^^Sit - t'), f scales as 
[f] The scaling dimensions of x and t can be determined 
by requiring the scale invariance of [d,x/f] = [x]/[f]'^^ and 
[J3x^/f] = [xf/(N[t]-^'^). The resulting scaling of [jc] = A^'^^ 
confirms our predictions based on the thermodynamic ap- 
proach (69). The scaling dimension of the time coordinate 
[f] - N^^^^ tells us that the spectrum is a scaling function of 
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N^^^ci). The present scaling analysis may also explain the re- 
cent numerical results of Ref. [49], where the eigen-frequency 
of the system at the critical point was found to scale as 1 /N"^ 
with e = 0.44 ± 0.02. It should be stressed that, remarkably, 
our scaling result differs from the one predicted for the fi- 
nite temperature transition of a closed system [56], for which 
[t] = N~^l'^ and the spectrum is a scaUng function of N^^'^oj. 



VI. KELDYSH APPROACH FOR ATOM OBSERVABLES 

We now analyze the single-atom observables of the open 
Dicke model. To this end, we represent each atom by a real 
field variable <p(, with ( = \,..,N. In order for <p( to be a 
faithful representation of an Ising spin, we constrain its fluc- 
tuations to lay on a unit Bloch sphere. This single-atom rep- 
resentation of the Ising spins allows to treat the qualitative 
effects of atom decay within a simplified "friction model" for 
(f)f which we explain below. This process couples directly to 
the local Ising degrees of freedom of the single atoms. The 
same is true for disorder due to spatial variations of the qubit- 
photon couplings [12]. In such cases, one cannot employ the 
single large-A^ spin representation of the Dicke model. 

We proceed by introducing into the path integral A' La- 
grange multipliers A(, corresponding to a suitable Fourier 
representation of the delta function [12, 50], 5(0^ - 1) = 

JcfAexp'''^'^^"''. On the closed time contour, this amounts to 
adding the following expression to the action: 



'-i,± 



\^'''^^^ ('^'•^^'^ " ^) " ^'--^^^ (^'--^'^ " ^) 



(70) 



Moving to the "classical/quantum" notation and adding the 
bare action of the atoms we obtain: 



(71) 



Note thai ^JA±[ can be associated with the excitation energy 
of the Ising spins and is to be determined self-consistently. 
The atom self-energies l^^^{a>) will be explained below. 
Finally, we have for the atom-cavity interaction: 

S<pa= y f K(0 + <») - 0-/(0 {a-{t) + a'_{t)) 1 . 

(72) 

The Keldysh action for the full Dicke model (1) then becomes 

S\_a,^,X\=Sa + S^^^x + ^^a , (73) 



with the various terms given by Eqs. (17,71, 72). 



We model local, single-atom damping in a simple effec- 
tive way which is consistent with symmetry properties of our 
real-valued Ising oscillators 0. The atoms are subject to decay 
into photon modes outside the cavity and possible other damp- 
ing mechanisms like s-wave scattering with other momentum- 
modes, trap loss or finite-size dephasing [51]. As a result 
some fraction of the atoms leave the two-density mode Hilbert 
space which maps to the Dicke model; others may be sponta- 
neously scattered back in. Representing the atoms by a com- 
plex field <5*, O, we subsume the above processes into Marko- 
vian decay of the atoms with the self-energies 



(74) 



with y an effective single-atom decay rate. Our effective 
real-valued Ising field in Eq. (71) may be viewed as the 
real component of the originally complex boson Oq/ci(0 = 

(0q/ci(O + J^q/ci(0)- Integrating out the ^-component, 

9 9 

A R K (xi~ + y + h)l 

ipi) = -ly^, ^4^) = +i7^> ^0 = ly — ' 

'(75) 



Note that this simple model for spontaneous emission couples 
to the cr^ projection of the atomic states and does not specify 
the states of the ctj projection of the spins. We emphasize, 
however, that the form of the dissipative self-energies is dic- 
tated by the combination of low-frequency expansion and the 
real-valued nature of the Ising field 0. In particular, a fre- 
quency independent term is ruled out for S^-^. The above re- 
sults for Markovian baths should be compared with the results 
for atoms in equilibrium. 



{co) = -iecj, 2^Eo(^) = 



2^£q(w) = 2ieo) coth 



,EQ' 

2r 



(76) 



where one also lets e — > at the end of the calculation. We 
here analyze the Dicke model in terms of the atomic degrees 
of freedom alone. One can exactly integrate out the pho- 
tons from the action (73). This is conveniently done by go- 
ing to a coordinate representation of the photons: aq/ci(0 = 

[xq/clit) + ipq/clit)), and first performing the integration 
over pq/ci and subsequently over Xq/d. We obtain, with S^^^^ 
given by Eq. (71), the atom-only action. 



S [<f>, A] - S 00,^ -1- S ^^^g. 

N ,2 



t,m=\ 



(77) 
(78) 



O-^(W) \ I 0cl,m(w) 



13 



where the matrix entries are 



-20Jn 



{oj + Ik) - oJq 

f.) 



/ ■ \2 71- 
(pJ-lK) - OJq\ 



(79) 



We now decouple Eq. (78) with a Hubbard-Stratonovich field 
il/e{io) <-> ^({(jS) and integrate out the ^-field. We assume tl/ to 



be time-independent, iA(w) = tff/(2n)6co,o, and substitute the 
Lagrange multipliers by a single effective field Aq/a/ — > Aq/d. 
The resulting Keldysh partition function 



DipDAe 



(80) 



obtains a prefactor of the number of atoms A'^ in the exponent 
multiplying the action 



(S^(0) - Aq) (cr'^(O))' + 2Aacr\0),/fq {cr\0)</fci + (r\0),J^q) - {o^iO)i/fci + cr\0)i/,qf 



4^(s^(0)-Aq) + 4 



Taking — > oo, we now extract the phase diagram, response 
and correlation functions, and the value of the order parame- 
ter using a saddle-point approximation. This can be obtained 
by requiring the derivatives with respect to Aq and tj/q to be 
zero, and then substituting Aq = 0, Ad = A, tfrq = 0, i^ci = 
The derivative with respect to Aq constrains -by construction 
Eq. (70)- the frequency integral of the Keldysh Green's func- 
tion to be equal to unity: 



with 



GL(a>) 



(82) 



2{oj^-A + T^ico)) {co^-A + ^(co)) 

2 2 

-2m^(w)^(crR(0))V. (83) 
The saddle-point condition for the order parameter yields 



dS 

dlffq 



= 



1 - 



tr«(0) 



= . 



(84) 



To determine the position of the Dicke transition, we need to 
compute the saddle point value of Ad = A m the normal and 
in the superradiant phases and equate the two values. In the 
normal (N) phase A is determined by (82) and (83) with \ff = Q, 



(85) 



Note that naively taking y ^ does not reproduce the equilib- 
rium value for Ajv, cf. also Subsec. VI C. In the ferromagnetic 
phase (FM) the order parameter acquires a finite expectation 
value \fi + ^ and, to fulfill Eq. (84), we need to require the 



argument of the square bracket to be zero: 



(86) 



where we have used Eq. (79) for cr'*(0). At the phase bound- 
aries both (85) and (86) must hold, leading to: 



8c 



(r2 + a>2)(/c2 + w2) 



3a)oWj. 



(87) 



Eqs. (82) and (84) determine the value of ferromagnetic order 
parameter x^j as well. When approaching the phase transition 
from above (g > gc). ^ vanishes as 



3(g+gc) r—— 
2 yg-gc 



\ 4^2 



(88) 



Compared with a closed system at equilibrium, Eq. (92), the 
order parameter for this open Markovian system vanishes with 
enhanced amphtude but with the same mean-field like square- 
root exponent. 



A. Atom spectral response function 

The single-atom retarded and advanced Green's function 
are determined by the derivative of (81) with respect to Ad 
and reads 



(89) 



2(6;2-A + S^(w)) ' 
from which follows the spectral response function 

Ji^^iu) = -2ImG^.(w) (90) 
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FIG. 6: (Color online) Atomic spectral response function Jl^^(a>) 
and atomic correlation function C^^icjS) in the normal phase (g < g^., 
green, dashed) and in the superradiant phase ig = i .2g,., black, solid). 
Other parameters used: cjq = co- = \, y = k = 0.2 (leading to 
gc ~ 0.6). The arrow illustrates the delta-function contributions from 
Eq. (83) in the superradiant phase. 



at w = 0. This is the non-equilibrium signature from sponta- 
neous emission of the pumped atoms into the vacuum outside 
the cavity (a continuum of modes with characteristic frequen- 
cies orders of magnitudes lower than the optical photons the 
pumped atoms emit when they spontaneously decay). 



C. Comparison with a closed system at equilibrium 

Note that taking k ^ and y — > in Eq. (87) does not 
reproduce the equilibrium value for a closed system. This 
is due to non-commuting limits of making the Markov ap- 
proximation and performing the integral to fulfill the sum rule 
Eq. (82). To obtain g^^ one needs to use the equilibrium bath 
self-energies Eqs. (28,76) from the start of the calculation and 
obtains the equilibrium analogs of Eqs. (87,88): 

1 



f Q = - 

^8 -8c 



(92) 



We note that, both dissipative channels, cavity photon loss and 
atomic spontaneous emission, shift the critical value of the 
coupling. The amplitude with which the feiTomagnetic order 
parameter vanishes is also different. 

Note that in a model for "one-way" spontaneous emission 
coupling to cr+ and cr" starting from a fully polarized atomic 
state as assumed in Sec. II A, one would recover the equi- 
librium limit for y ^ 0. As explained above, our "mixing" 
model for spontaneous emission couples to cr' and thereby 
assumes a mixed state of the atoms (similar to a many-body 
paramagnet). 



D. Atom distribution function and low-frequency effective 
temperature 



as the expected frequency-resolved signal from the atoms after 
local, time-modulated density perturbations. Fig. 6 displays 
the characteristic Lorentzian shape of the spectral response 
function peaked at V^, broadened by single atom decay ~ 
y. The single atom response is smooth across the transition. 
This is not to be confused with the roton-type mode softening 
observed by Mottl et al. [51] that pertains to the collective 
atomic density excitations for a finite number of atoms. 



We now execute the procedure of subsection HID to cal- 
culate the effective temperature of the atoms. With the atom 
Green's functions and the simple model for atom decay pre- 
sented above, one finds 



+ + dLt? 1 



2a), 



(93) 



B. Atom correlation function 

The atom correlation function is given by the Keldysh 
Greens function Eq. (83) 



(91) 



G^^ is defined in (83), and is exhibited in Fig. 6. Even before 
the onset of the superradiance peaks (black and blue-dashed 
arrows) for g > gc, the correlation function has finite weight 



leading to the effective temperature 



neff. r'+W? 



4<jj. 



(94) 



which is independent of the coupling strength to the photons 
(cf. also Fig. 7). For y <K cj,, the effective temperature is 
set by the recoil energy of the atoms Er - co^jl. Within our 
model, T'^ also does not depend on the cavity loss rate k, con- 
trary to what obtains for the model of Ref. 13. In our case the 
reason for this is the careful treatment of the thermodynamic 
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FIG. 7: (Color online) Distribution function F^^iuj) for the same nu- 
merical parameters as Fig. 6. The expression for the distribution 
function Eq. (94) is independent of the atom-photon coupling g. 

limit N oo limit leading to Eqs. (77,81). This limit ensures 
that the only photon-induced self-energies for the atoms occur 
for zero-frequency quantities (the weight of 5((jS) in Eq. (83)). 



VII. CONCLUSION 

In this paper, we presented a path integral approach for the 
non-equilibrium steady-states of driven quantum systems cou- 
pled to Markovian baths, such as ultracold atoms in optical 
cavities. In the past, these systems have more often been de- 
scribed using a Master equation formalism. We believe that 
our Keldysh approach allows an easier comparison with other 
equilibrium and non-equilibrium (classical and quantum) sys- 
tems. While some of the results presented here are actually 
new, and not just known results re-phrased in a new approach, 
the full utility of our approach will become clear then com- 
puting thermodynamics and critical properties of large, open 
systems with spatially fluctuating degrees of freedom such as 
disorder [12, 14]. In these correlated quantum many-body 
situations. Master equation approaches are typically limited 
to relatively small number of atoms and a recipe to compute 
disorder-averaged quantities does not seem to exist. 

We first applied our formalism to the cavity vacuum and 
subsequently added atomic qubits, interacting with the cavity 
through a Dicke interaction, and computed the key observ- 
ables for both, the photons and atoms. In the presence of 
the atoms, the cavity exhibits a thermodynamic "superradi- 
ant" phase transition, whose critical properties we analyzed in 
detail. They differ from the zero-temperature quantum phase 
transition of the corresponding closed system due to the pres- 
ence of an emergent finite "effective temperature" for the rel- 
evant low-energy photons, induced by the atom-photon inter- 
action. Due to the non-equilibrium nature of the problem, 
the effective temperature characterizes only the low-frequency 
behavior of the system. At higher frequencies, the distribu- 
tion functions display non-equilibrium and quantum behav- 
iors. The definition of an effective temperature is observable 



dependent, leading to different effective temperatures for the 
atoms versus the photons. 

In the future, it will be interesting to apply our approach 
to dissipative quantum glasses coupled to Markovian (and 
other) baths such as potentially achievable in multi-mode op- 
tical cavities [12-14] or circuit QED [52]. It would also be 
desirable to obtain a more general classification of conditions 
under which quantum phase transitions of closed systems are 
turned into thermal phase transitions by dissipation-and per- 
haps to find counterexamples by engineered dissipation along 
the lines ofRefs. [38, 39]. 
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Appendix A: Self-energy of an open cavity in tlie rotating frame 

In this appendix, we discuss how the coherent drive with a 
frequency scale Up, which exceeds all other frequency scales, 
justifies the form of the Markovian dissipative action Eq. (16), 
which in particular displays frequency independent terms only 
which are 5-coiTelated in time and neglect memory effects. 

Our starting point is the Hamiltonian of a single boson a, 
coupled to a continuum of vacuum fields i^^, via 

Hq = a)ca'^a + ^a)kil/[>l/k +gk(a}ipk + a>lj[)- (Al) 
k 

Here gk is the coupling constant between the cavity boson and 
the external vacuum, and we neglected counter-propagating 
terms of the form aV^- Eq. (Al) is quadratic in i/r^, allowing 
us to analytically integrate-out the vacuum fields and obtain 
a cavity-only action of the form (17). If we assume that the 
vacuum fields are kept at an equilibrium temperature Text = 
300/r, the corresponding entries are 

[G*(w)]"' = u-u)c-5u) + iK{u>), 

= 2//:(w)cothJ^j. (A2) 

Here 6u) corresponds to the Lamb shift and can be absorbed 
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in a finite renormalization of ojc- The function Kioj) = 
Tik \8l\^i^ - (^k) is the spectral density of the vacuum. 

The inverse Green's functions (A2) describe the cavity 
mode in the lab frame. In practice, it is often more convenient 
to move to a frame rotating with a constant frequency, in our 
case corresponding to the pump frequency ojp. In this frame, 
the photons are described by Eq. (A2) with cj cOp + (o. 
Next, we apply the equivalent of the Wigner-Weisskopf ap- 
proximation, by Taylor expanding the inverse Green's func- 
tion in small o) to zero order. This approximation is justified 
by the energy scale separation discussed in the text. We obtain 

[G*(a>)]~^ = CO - ojc - Sco - ojp + iK{a>p), 

D*^ = 2iK{up)coihi^^^. (A3) 

The factor coth^^^) plays the role of the 2n + 1 factor ap- 
pearing in the finite temperature extension of the master equa- 
tion (2). To be precise, the two expressions coincide only for 
a>c = (Op- For most experiments the pump frequency is any- 
way much higher than the external temperature and we can 
approximate coth^j^) = 1- Under this approximation, Eq. 
(A3) becomes equivalent to Eqs. (20-21) with 



oJr+Soj- a>„ 



K{cOp). 



(A4) 



Appendix B: Photon-only action in tlie superradiant pliase 



In the superradiant (SR) phase, the field a (b) is given 
by the sum of a time-independent component ao (bo) and a 
fluctuating term. The action governing the fluctuating terms 
can be obtained from Eq. (35) by substituting a ^ ao + 6a 
(b bo+ 6b) (we choose uq, bo real without loss of general- 
ity). At ^ 00, we end up with the quadratic action 



>SR 





[G4X4] 



(CO) 



[G4X4] 

A- 



6yi(ui) 



(Bl) 



with the Green's functions [^4^4] '(w) - ([^4x4] '(w)^ 



a>o -1- IK 





-g 


-g 





-0) -0)0- IK 


-g 


-g 


-g 


-g 


oj- 0}^- 26(0^ 


-6cOz 


-g 


-g 


-6a>, 


-U)-0}^- 



Here the 8 -vector SViipi) is defined in the analogous way to 
Y^ioS) ofEq. (40) and 



3^ 2 3^2+ 2 
g = g- -ttK ~ 



4g 



(B3) 



S'iOz = ^boiao + a;) » -4aj,^--^. (B4) 

We note that the principal change to the spectral response 
and correlation fimction in the superradiant phase 

Jlaa^ioj) = JlgaM^) , (B5) 
Caa^Oj) = CsaMi^) + ^0^oj,0 (B6) 

is the 5-function peak at w = in the correlation function due 
to coherent photons ("photon condensate"). 



Appendix C: Damped dynamics near tlie pliase transition 

We argue based on a systematic low-frequency expansion 
of the inverse retarded Green's function that the overdamped 
dynamics observed in the vicinity of the phase transition is 
generic for systems where a phase transition is driven by a 
competition within the Hamiltonian sector, while dissipative 
dynamics acts as a "spectator". To see this, we (i) write the 
most general form of the inverse retarded Green's function 



pi(o) 0i(0) 
o-'i-co) p'i-o>) 



)• 



(CI) 



and (ii) use that the phase transition is governed by low fre- 
quency behavior and an expansion in powers of the frequency 

is appropriate, 

p(oj) = -V + zoJ, o(oj) = -ju + yoj (C2) 
with complex coefficients and a low frequency spectrum 
ilm[z*y - y'li] ± V(kP - \yPW - (Im[z*v - y*fi])^ 



0J± 



(kP - lyP) 



(C3) 



Without dynamic renormaUzation effects, z = 1 and 
all other frequency coefficients are zero, so they will 
be generically much smaller than one (more precisely, 
|Im[z]|,|Re[>']|,|Im[>']| «: 1), and in particular |zp » Ijp. (In 
the large A'^ open Dicke model, they are exactly zero.) A mass 
gap, i.e. the scale that characterizes the action at zero fre- 
quency, provides a measiu'e of the distance from the phase 
transition and reads 



= detGf-2 = 0) = \vf - > 



(C4) 



= 2i diag(/c, k, 0, 0). 



(B2) 



(the last inequahty must hold for a stable physical system). 
Approaching the phase transition, this gap shrinks to zero, 
such that the frequencies must become purely imaginary as 
a generic feature of a phase transition in the presence of dis- 
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sipation. Indeed, in a situation where the phase transition is 
driven by a competition within the Hamiltonian sector of the 
problem by a quantity g, the dominant g dependence is con- 
tained in a^ig) (more precisely, Re[v],Re|ju]), while the dissi- 
pative scales (Im[v],lm[yuJ) do not strongly depend on g and 
remain essentially at their bare, finite values even at the tran- 
sition point. (In the open Dicke model, only the real parts are 
modified, while the imaginary parts exactly remain at their 
bare values.) In such a situation, as a{g — > gc) — > 0, we may 
expand the square root in Eq. (C3) in a^, identifying that pa- 
rameter as the distance from the phase transition. 

Appendix D: divergence in Marlsov distribution functions 

We here show that the l/oj pole in the photon distribution 
function at low frequency, and the associated low-frequency 
efi'ective temperature (LET), is indeed a generic feature of 
Markovian non-equilibrium systems. 

To this end, we consider the low-frequency regime, where 
in the spirit of a systematic derivative expansion the inverse 
retarded Green's function of the photon takes the form 

G*"' =(wH-iV2)crj-//, H ^ Vil + IJiO-jt + IU2Cry, (Dl) 



where vi,//i(y2,yU2) denote the real (imaginary) part of v,/!. 
We set z = l,y = here, which in principle contribute at 
0{(jl)), and anticipate that this omission will not alter the qual- 
itative results. The hermitean part H represents Hamiltonian 
dynamics, the antihermitean ~ icr^ decay. In a derivative ex- 
pansion, the most general form of the Keldysh component is 

= 2i(/fil -I- K2crjc). (D2) 

Solving the fluctuation-dissipation relation, we obtain 

Ficj) = —a^ - -\ — {ni(Tx+H2(Ty) + -I- K2(tX (D3) 



Crucially, this confirms the Ijcj divergence behavior of the 
distribution fimction. Allowing for the most general form of 
of the frequency expansion in terms of a finite imaginary 
part of z and finite y only results in subleading corrections: 
also in this case, lim„^o['JJ ■ Fid))] const. Clearly, adding 
frequency dependent terms to D'^ only leads to subleading 
corrections in F. Therefore, the Ijco pole at low frequency, 
and the associated scale generated in this regime, the LET, is 
a generic feature of Markovian non-equilibrium systems. 
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